function f = z_high(zh,beta,k,R,r0,q_bar)
% Upper bound of branches leading to bunching
f = beta*exp((R - beta - r0)/beta)*zh - k -q_bar*(R - r0 - beta*log(q_bar) + beta*log(zh));